home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Cream of the Crop 21
/
Cream of the Crop 21 (Terry Blount) (October 1996).iso
/
os2
/
e33el2.zip
/
emacs
/
19.33
/
lisp
/
float.el
< prev
next >
Wrap
Lisp/Scheme
|
1996-01-20
|
16KB
|
459 lines
;;; float.el --- floating point arithmetic package.
;; Copyright (C) 1986 Free Software Foundation, Inc.
;; Author: Bill Rosenblatt
;; Maintainer: FSF
;; Keywords: extensions
;; This file is part of GNU Emacs.
;; GNU Emacs is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 2, or (at your option)
;; any later version.
;; GNU Emacs is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
;; GNU General Public License for more details.
;; You should have received a copy of the GNU General Public License
;; along with GNU Emacs; see the file COPYING. If not, write to the
;; Free Software Foundation, Inc., 59 Temple Place - Suite 330,
;; Boston, MA 02111-1307, USA.
;;; Commentary:
;; Floating point numbers are represented by dot-pairs (mant . exp)
;; where mant is the 24-bit signed integral mantissa and exp is the
;; base 2 exponent.
;;
;; Emacs LISP supports a 24-bit signed integer data type, which has a
;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
;; This gives six significant decimal digit accuracy. Exponents can
;; be anything in the range -(2**23) to +(2**23)-1.
;;
;; User interface:
;; function f converts from integer to floating point
;; function string-to-float converts from string to floating point
;; function fint converts a floating point to integer (with truncation)
;; function float-to-string converts from floating point to string
;;
;; Caveats:
;; - Exponents outside of the range of +/-100 or so will cause certain
;; functions (especially conversion routines) to take forever.
;; - Very little checking is done for fixed point overflow/underflow.
;; - No checking is done for over/underflow of the exponent
;; (hardly necessary when exponent can be 2**23).
;;
;;
;; Bill Rosenblatt
;; June 20, 1986
;;
;;; Code:
;; fundamental implementation constants
(defconst exp-base 2
"Base of exponent in this floating point representation.")
(defconst mantissa-bits 24
"Number of significant bits in this floating point representation.")
(defconst decimal-digits 6
"Number of decimal digits expected to be accurate.")
(defconst expt-digits 2
"Maximum permitted digits in a scientific notation exponent.")
;; other constants
(defconst maxbit (1- mantissa-bits)
"Number of highest bit")
(defconst mantissa-maxval (1- (ash 1 maxbit))
"Maximum permissible value of mantissa")
(defconst mantissa-minval (ash 1 maxbit)
"Minimum permissible value of mantissa")
(defconst floating-point-regexp
"^[ \t]*\\(-?\\)\\([0-9]*\\)\
\\(\\.\\([0-9]*\\)\\|\\)\
\\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
"Regular expression to match floating point numbers. Extract matches:
1 - minus sign
2 - integer part
4 - fractional part
8 - minus sign for power of ten
9 - power of ten
")
(defconst high-bit-mask (ash 1 maxbit)
"Masks all bits except the high-order (sign) bit.")
(defconst second-bit-mask (ash 1 (1- maxbit))
"Masks all bits except the highest-order magnitude bit")
;; various useful floating point constants
(setq _f0 '(0 . 1))
(setq _f1/2 '(4194304 . -23))
(setq _f1 '(4194304 . -22))
(setq _f10 '(5242880 . -19))
;; support for decimal conversion routines
(setq powers-of-10 (make-vector (1+ decimal-digits) _f1))
(aset powers-of-10 1 _f10)
(aset powers-of-10 2 '(6553600 . -16))
(aset powers-of-10 3 '(8192000 . -13))
(aset powers-of-10 4 '(5120000 . -9))
(aset powers-of-10 5 '(6400000 . -6))
(aset powers-of-10 6 '(8000000 . -3))
(setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
highest-power-of-10 (aref powers-of-10 decimal-digits))
(defun fashl (fnum) ; floating-point arithmetic shift left
(cons (ash (car fnum) 1) (1- (cdr fnum))))
(defun fashr (fnum) ; floating point arithmetic shift right
(cons (ash (car fnum) -1) (1+ (cdr fnum))))
(defun normalize (fnum)
(if (> (car fnum) 0) ; make sure next-to-highest bit is set
(while (zerop (logand (car fnum) second-bit-mask))
(setq fnum (fashl fnum)))
(if (< (car fnum) 0) ; make sure highest bit is set
(while (zerop (logand (car fnum) high-bit-mask))
(setq fnum (fashl fnum)))
(setq fnum _f0))) ; "standard 0"
fnum)
(defun abs (n) ; integer absolute value
(if (>= n 0) n (- n)))
(defun fabs (fnum) ; re-normalize after taking abs value
(normalize (cons (abs (car fnum)) (cdr fnum))))
(defun xor (a b) ; logical exclusive or
(and (or a b) (not (and a b))))
(defun same-sign (a b) ; two f-p numbers have same sign?
(not (xor (natnump (car a)) (natnump (car b)))))
(defun extract-match (str i) ; used after string-match
(condition-case ()
(substring str (match-beginning i) (match-end i))
(error "")))
;; support for the multiplication function
(setq halfword-bits (/ mantissa-bits 2) ; bits in a halfword
masklo (1- (ash 1 halfword-bits)) ; isolate the lower halfword
maskhi (lognot masklo) ; isolate the upper halfword
round-limit (ash 1 (/ halfword-bits 2)))
(defun hihalf (n) ; return high halfword, shifted down
(ash (logand n maskhi) (- halfword-bits)))
(defun lohalf (n) ; return low halfword
(logand n masklo))
;; Visible functions
;; Arithmetic functions
(defun f+ (a1 a2)
"Returns the sum of two floating point numbers."
(let ((f1 (fmax a1 a2))
(f2 (fmin a1 a2)))
(if (same-sign a1 a2)
(setq f1 (fashr f1) ; shift right to avoid overflow
f2 (fashr f2)))
(normalize
(cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
(cdr f1)))))
(defun f- (a1 &optional a2) ; unary or binary minus
"Returns the difference of two floating point numbers."
(if a2
(f+ a1 (f- a2))
(normalize (cons (- (car a1)) (cdr a1)))))
(defun f* (a1 a2) ; multiply in halfword chunks
"Returns the product of two floating point numbers."
(let* ((i1 (car (fabs a1)))
(i2 (car (fabs a2)))
(sign (not (same-sign a1 a2)))
(prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
(lohalf (* (hihalf i1) (lohalf i2)))
(lohalf (* (lohalf i1) (hihalf i2)))))
(prodhi (+ (* (hihalf i1) (hihalf i2))
(hihalf (* (hihalf i1) (lohalf i2)))
(hihalf (* (lohalf i1) (hihalf i2)))
(hihalf prodlo))))
(if (> (lohalf prodlo) round-limit)
(setq prodhi (1+ prodhi))) ; round off truncated bits
(normalize
(cons (if sign (- prodhi) prodhi)
(+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
(defun f/ (a1 a2) ; SLOW subtract-and-shift algorithm
"Returns the quotient of two floating point numbers."
(if (zerop (car a2)) ; if divide by 0
(signal 'arith-error (list "attempt to divide by zero" a1 a2))
(let ((bits (1- maxbit))
(quotient 0)
(dividend (car (fabs a1)))
(divisor (car (fabs a2)))
(sign (not (same-sign a1 a2))))
(while (natnump bits)
(if (< (- dividend divisor) 0)
(setq quotient (ash quotient 1))
(setq quotient (1+ (ash quotient 1))
dividend (- dividend divisor)))
(setq dividend (ash dividend 1)
bits (1- bits)))
(normalize
(cons (if sign (- quotient) quotient)
(- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
(defun f% (a1 a2)
"Returns the remainder of first floating point number divided by second."
(f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
;; Comparison functions
(defun f= (a1 a2)
"Returns t if two floating point numbers are equal, nil otherwise."
(equal a1 a2))
(defun f> (a1 a2)
"Returns t if first floating point number is greater than second,
nil otherwise."
(cond ((and (natnump (car a1)) (< (car a2) 0))
t) ; a1 nonnegative, a2 negative
((and (> (car a1) 0) (<= (car a2) 0))
t) ; a1 positive, a2 nonpositive
((and (<= (car a1) 0) (natnump (car a2)))
nil) ; a1 nonpos, a2 non